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Executive  Summary 


The  Weather  Research  and  Forecasting  Model  (WRF)  is  a  numerical  weather 
prediction  model  that  has  been  used  for  many  applications,  including  My  Weather 
Impacts  Decision  Aid  (MyWIDA).  WRF  is  maintained  by  the  National  Center  for 
Atmospheric  Research  (NCAR),  which  has  developed  a  suite  of  Model 
Evaluation  Tools  (MET)  to  evaluate  the  accuracy  of  WRF  forecasts  using 
observations  of  meteorological  variables  such  as  temperature,  relative  humidity, 
and  wind.  The  traditional  use  of  MET  calculates  model  performance  over  the 
entire  model  domain. 

For  Soldiers  in  theater,  the  Army  requires  high-resolution  weather  forecasting  to 
resolve  atmospheric  features  with  wavelengths  on  the  order  of  5  kilometers  (km) 
or  less.  This  requires  models  that  operate  on  a  model  grid  spacing  of  about  1  km 
or  less  in  the  finest,  or  most  resolved,  domain.  To  validate  such  high-resolution 
modeling  against  observations  requires  a  more  focused  spatial  and  temporal 
approach  over  parts  of  the  domain  rather  than  the  domain  as  a  whole.  With  a 
Geographic  Information  System  (GIS),  researchers  can  now  consider  terrain 
type/slope,  land  use  effects,  and  other  spatial  and  temporal  variables  as 
explanatory  metrics  in  model  assessments.  GIS  techniques,  when  coupled  with 
high-resolution  point  and  gridded  observation  sets,  allow  location-based 
approaches  that  pennit  discovery  of  spatial  and  temporal  scales  where  models  do 
not  sufficiently  resolve  the  desired  phenomena — for  example,  turbulence  effects 
or  mountain  and  lee  waves. 

In  this  technical  report,  we  discuss  our  initial  efforts  in  using  GIS  tools  to  study 
model  errors  in  the  Advanced  Research  WRF  (WRF-ARW)  with  a  1  -km 
horizontal  grid  spacing  inner  domain  centered  near  San  Diego,  California.  The 
San  Diego  area  contains  a  mixture  of  urban,  suburban,  agricultural,  and 
mountainous  terrain  types  along  with  a  rich  array  of  observational  data  with 
which  to  illustrate  our  ability  to  conduct  subdomain  verification.  We  selected  5 
case  study  days  from  February  and  March  of  2012  with  varied  synoptic 
conditions. 

A  literature  review  indicated  that  elevation  is  strongly  correlated  with 
meteorological  parameters,  at  least  in  a  climatic  sense;  consequently  we  focused 
this  report  on  elevation  as  the  explanatory  variable.  We  found  that  elevation 
accounts  for  a  significant  portion  of  the  variance  in  the  model  error  of  surface 
temperature  and  relative  humidity  predictions.  On  average,  elevation  accounts  for 
about  10%  to  20%  of  the  total  variance.  In  some  cases  however,  elevation 
explained  more  than  50%  of  the  total  variance,  while  in  other  cases  elevation  was 
not  a  significant  explanatory  variable  at  all.  Overall  we  found  that  elevation 
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accounts  for  more  variance  during  the  afternoon  hours  and  less  variance  during 
the  early  morning  hours;  however,  we  could  not  discern  a  clear  diurnal  pattern 
within  the  model  errors. 

In  this  study,  we  demonstrated  the  effectiveness  of  using  a  GIS  to  evaluate  model 
perfonnance  at  the  subdomain  level.  A  GIS  can  analyze  spatially  distributed  point 
forecast  errors  and  show  the  dependence  of  forecast  errors  on  terrain 
characteristics  at  the  spatial  scales  of  phenomena  of  interest  to  Warfighters,  such 
as  mountain/valley  breezes  and  land/sea  breezes.  The  GIS  techniques  used  in  this 
study  show  considerable  promise  for  demonstrating  the  relationships  of  high- 
resolution  model  errors  to  explanatory  variables  such  as  terrain  as  well  as 
identifying  subdomains  where  the  spatial  statistics  are  more  homogeneous. 
Further  analysis  using  other  explanatory  variables  such  as  slope,  land  use,  and 
distance  from  the  coastline  will  be  needed  to  draw  firm  conclusions  about  the 
model  performance.  We  achieved  our  goal  with  this  study,  which  was  to  develop 
a  GIS  analysis  method  that  goes  beyond  the  types  of  analysis  and  conclusions 
gained  from  domain-level,  aggregate  statistics. 


1.  Background 


Weather  has  a  significant  impact  on  Army  personnel,  weapons,  tactics,  and 
operations;  therefore,  accurate  weather  forecasts  can  be  a  deciding  factor  in  any 
conflict — large  or  small.  A  famous  example  is  D-day  in  World  War  II  when  the 
Allied  nations  invaded  Normandy,  France.  Only  a  few  days  each  month  had  the 
timing  of  high  and  low  tides  and  full  moon  conditions  suitable  for  landing 
thousands  of  men  on  the  beaches  of  Normandy.  On  June  4,  1944,  high  winds  and 
heavy  seas  made  a  landing  impossible;  however,  James  Stagg  of  the  Royal  Air 
Force  and  his  team  of  meteorologists  predicted  that  the  weather  would  improve  on 
June  5,  and  General  Eisenhower  decided  to  go  ahead  with  the  invasion.  The 
German  meteorologists  were  not  expecting  the  weather  to  break  and  were  thus 
surprised  by  the  invasion  (Ross  2014). 

As  computing  technology  has  advanced,  the  weather  forecasting  task,  once  the 
primary  role  of  a  human  forecaster  in  theater,  has  shifted  to  computerized 
Numerical  Weather  Prediction  (NWP)  models.  Scientists  around  the  world  have 
used  the  Weather  Research  and  Forecasting  (WRF)  model  extensively  for  many 
applications.  In  this  study,  we  have  used  the  Advanced  Research  version  of  WRF 
(Skamarock  et  al.  2008)  that  we  abbreviate  as  WRF-ARW.  WRF-ARW  includes 
Four  Dimensional  Data  Assimilation  (FDDA)  techniques  that  can  be  used  to 
incorporate  observations  into  the  model  so  that  forecast  quality  is  improved  (Deng 
et  al.  2009;  Stauffer  and  Seaman  1994).  The  US  Army  Research  Laboratory 
(ARL)  uses  WRF-ARW  as  the  core  of  its  Weather  Running  Estimate-Nowcast 
(WRE-N)  weather  forecasting  model. 

The  Army  requires  high-resolution  weather  forecasting  to  resolve  atmospheric 
features  with  wavelengths  on  the  order  of  5  km  or  less,  which  imposes  a 
requirement  for  NWP  to  operate  on  a  model  grid  spacing  on  the  order  of  1  km  or 
less  in  the  finest,  or  most  resolved,  domain  to  resolve  weather  phenomena  of 
interest  to  the  Soldier  in  theater.  The  atmospheric  flows  of  interest  to  the  Army 
include  mountain/valley  breezes,  sea  breezes,  and  other  flows  induced  by 
differences  in  land  surface  characteristics.  High-resolution  NWP  forecasts  need  to 
be  validated  against  observations  before  their  outputs  can  be  used  by  applications 
such  as  My  Weather  Impacts  Decision  Aid  (MyWIDA)  developed  by  Brandt  et  al. 
(2013).  Weather  forecast  validation  has  always  been  of  interest  to  the  civilian  and 
military  weather  forecasting  community;  see  for  example  the  reviews  by  Ebert  et 
al.  (2013)  and  Casati  et  al.  (2008),  or  the  guides  by  Jolliffe  and  Stephenson  (2011) 
or  Wilks  (2011).  The  validation  of  the  models,  especially  high-resolution  NWP, 
has  proven  to  be  especially  difficult  when  addressing  small  temporal  and  spatial 
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scales  (National  Research  Council  2010)  that  characterize  NWP  for  use  in  Army 
applications. 

The  WRF  model  is  maintained  by  the  National  Center  for  Atmospheric  Research 
(NCAR),  which  has  also  developed  a  suite  of  Model  Evaluation  Tools  (MET) 
(NCAR  2013)  to  evaluate  WRF-ARW  performance.  MET  Point-Stat  perfonns 
traditional  grid-to-point  verification,  while  MET  Grid-Stat  performs  grid-to-grid 
neighborhood  verification  to  account  for  the  uncertainty  inherent  in  high- 
resolution  forecasting.  The  MET  Method  for  Object-based  Diagnostic  Evaluation 
(MODE)  has  been  used  to  develop  techniques  for  object-based  spatial  verification 
of  high-resolution  forecast  grids  of  continuous  meteorological  variables. 

ARL  has  employed  MET  in  prior  assessments  such  as  that  of  Raby  et  al.  (2012) 
who  evaluated  2  models  to  arrive  at  domain  level  conclusions  about  the  various 
strengths  and  weaknesses  of  these  models  and  their  accuracies.  While  MET 
proved  useful  as  an  assessment  tool  for  forecasts  over  a  regional  domain,  they 
found  its  ability  to  illuminate  model  shortcomings  on  the  spatial  and  temporal 
scales  needed  by  battlefield  commanders  to  be  lacking.  In  another  assessment 
Dumais  et  al.  (2012)  evaluated  3  models  that  produced  3-dimensional  (3-D)  wind 
fields  to  arrive  at  domain-level  conclusions.  Although  Dumais  et  al.  did  not 
employ  any  of  the  MET  in  their  analysis,  their  desire  was  to  also  identify  model 
strengths  and  weaknesses,  and  again,  they  found  their  ability  to  draw  spatial  and 
temporal  conclusions  lacking. 

In  other  related  work,  Johnson  (2012)  conducted  a  location-based  analysis  of 
WRF  over  an  inland  testing  range  and  for  specific  model  time  bins  for  a  single 
day  of  interest.  In  his  analysis,  Johnson  categorized  the  terrain  as  belonging  to 
Valleys,  Plains,  or  Mountains  and  concluded  that  this  categorization  increased 
correlation  coefficients  for  these  categories.  However,  the  differences  when 
considering  horizontal  grid  resolutions  of  1  km  versus  3  km  were  present  only  in 
the  “Valley”  and  “Mountain”  cases.  Though  differences  in  1  km  and  3  km  were 
present  (favoring  1  km  improved  accuracy),  they  were  minimal.  Johnson  et  al. 
(2014)  extends  and  refines  this  result. 

Our  assertion  (Smith  et  al.  2014a;  Smith  et  al.  2014b;  2015;  Smith  et  al.  2014c; 
Smith  et  al.  2013)  is  that  high-resolution  models  require  verification  on  temporal 
and  spatial  scales  appropriate  for  the  phenomena  that  are  being  forecast.  Our 
principal  thesis  is  that  a  Geographic  Infonnation  System  (GIS)  has  tools  that 
allow  us  to  divide  the  domain  into  more  homogenous  subdomains  according  to 
explanatory  variables  such  as  terrain  and  land  use.  Our  ultimate  goal  is  to  make 
more  accurate  inferences  about  model  performance  as  well  as  providing  a  better 
understanding  of  model  strengths  and  weaknesses. 
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The  purpose  of  this  study  is  to  develop  GIS  tools  that  can  better  evaluate  model 
performance.  Our  study  domain  is  located  in  the  San  Diego,  California  region, 
which  we  chose  because  of  its  varied  terrain  and  large  number  of  weather 
observing  stations.  We  chose  5  case  study  days  with  varied  synoptic  conditions 
from  February  and  March  of  2012.  We  used  the  MET  Point-Stat  to  calculate  the 
model  forecast  value  at  each  of  the  observing  stations  and  analyzed  the  matched 
forecast-observation  pairs  using  GIS  tools. 

2.  Domain  and  Model 


The  ARL  WRE-N  (Dumais  et  al.  2013;  Dumais  et  al.  2004)  has  been  designed  as 
a  convection-allowing  application  of  the  WRF-ARW  model  (Skamarock  et  al. 
2008)  with  an  observation-nudging  FDDA  option  (Deng  et  al.  2009;  Liu  et  al. 
2005).  For  this  investigation,  we  configured  WRE-N  in  a  multinest  (i.e.,  9/3/1 
km)  configuration  to  produce  a  fine  inner  mesh  with  1-km  grid  spacing  and 
leveraged  an  external  global  model  for  cold-start  initial  conditions  and  time- 
dependent  lateral  boundary  conditions  for  the  outermost  nest.  For  ARL 
development  and  testing,  this  global  model  has  been  the  National  Center  for 
Environmental  Prediction’s  Global  Forecast  System  (GFS)  model  (Environmental 
Modeling  Center  2003).  The  WRE-N  is  envisioned  to  be  a  rapid-update  cycling 
application  of  WRF-ARW  with  FDDA  and  optimally  could  refresh  itself  at 
intervals  up  to  hourly  (dependent  upon  the  observation  network)  (Dumais  et  al. 
2012;  Dumais  and  Reen  2013). 

For  this  particular  experiment,  the  model  was  run  with  a  base  time  of  1200 
coordinated  universal  time  (UTC)  and  generated  output  for  each  hour  from  1200 
UTC  to  1200  UTC  of  the  following  day  for  a  total  of  25  model  outputs  on  each  of 
5  days  in  February  and  March  of  2012.  The  modeling  domain  is  depicted  in 
Fig.  1. 
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Fig.  1  Triple-nested  model  configuration 

Note:  The  outer  domain  is  dimensioned  as  1,566  km  by  1,566  km  at  a  grid  spacing  of  9  km;  a 
middle  nest  domain  dimensioned  as  720  km  by  720  km  at  a  grid  spacing  of  3  km;  and  an 
inner  domain  dimensioned  as  126  km  by  126  km  at  a  grid  spacing  of  1  km.  Each  domain 
center  is  coincident  and  the  nested  configuration  is  centered  near  San  Diego,  California. 

2.1  Observations 

The  initial  conditions  were  constructed  by  starting  with  the  GFS  data  as  the  first 
guess  for  an  analysis  using  observations.  Most  observations  were  obtained  from 
the  Meteorological  Assimilation  Data  Ingest  System  (MADIS)  (National  Oceanic 
and  Atmospheric  Administration  [NOAA]  2014a),  except  for  the  Tropospheric 
Airborne  Meteorological  Data  Reporting  (TAMDAR)  (Daniels  et  al.  2004) 
observations,  which  were  obtained  from  AirDat  (2014).  The  MADIS  database 
included  standard  surface  observations,  mesonet  surface  observations,  maritime 
surface  observations,  wind  profiler  measurements,  rawinsonde  soundings,  and 
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Aircraft  Communications  Addressing  and  Reporting  System  (ACARS)  data.  Use 
and  reject  lists  were  obtained  from  developers  of  the  Real-Time  Mesoscale 
Analysis  system  (NOAA  2014b),  and  these  were  used  to  filter  MADIS  mesonet 
observations.  This  quality  assurance  evaluation  is  especially  important  given  the 
greater  tendency  of  mesonet  observations  to  be  poorly  sited  than  other,  more 
standard,  surface  observations. 

The  Obsgrid  component  of  WRF  was  used  for  quality  control  of  all  observations. 
This  included  gross  error  checks,  comparing  observations  to  a  background  field 
(here  GFS),  and  comparing  observations  to  nearby  observations.  We  modified 
Obsgrid  to  allow  for  single-level  observations  such  as  the  TAMDAR  and  ACARS 
data  to  be  more  effectively  compared  against  the  GFS  background  field.  We 
employed  observation  nudging  to  the  observations  from  these  same  sources  from 
12  to  18  UTC,  followed  by  1-h  ramping  down  of  the  nudging  from  18  to  19  UTC 
during  which  time  no  new  observations  are  assimilated.  The  forecast  period  thus 
begins  at  18  UTC  because  no  observations  after  this  time  are  assimilated. 

2.2  Parameterizations 

For  the  parameterization  of  turbulence  in  WRE-N,  a  modified  version  of  the 
Mellor-Yamada-Janjic  (MYJ)  Planetary  Boundary  Layer  (Janjic  1994)  scheme 
was  used.  This  modification  decreases  the  background  turbulent  kinetic  energy 
and  alters  the  diagnosis  of  the  boundary  layer  depth  used  for  model  output  and 
data  assimilation  (Reen  et  al.  2014).  The  WRF  single-moment,  5-class 
microphysics  parameterization  is  used  on  all  domains  (Hong  et  al.  2004),  while 
the  Kain-Fritsch  (Kain  2004)  cumulus  parameterization  is  used  only  on  the  9-km 
outer  domain.  For  radiation,  the  Rapid  Radiative  Transfer  Model  (RRTM) 
parameterization  (Mlawer  et  al.  1997)  is  used  for  long  wave  radiation  and  the 
Dudhia  (1989)  scheme  for  shortwave  radiation.  The  Noah  land  surface  model 
(Chen  and  Dudhia  2001a;  2001b)  is  used.  Additional  references  and  other  details 
for  these  parameterization  schemes  are  available  from  Skamarock  et  al.  (2008). 

2.3  Case  Study  Days 

The  case  study  days  were  selected  on  the  basis  of  the  prevailing  synoptic  weather 
conditions  over  the  nested  domains,  and  a  short  description  of  these  conditions  is 
provided  in  Table  1. 
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Table  1  Synoptic  conditions  for  the  case  study  days  considered 


Case 

Dates  (all  2012) 

Description 

1 

February  07-08 

Upper-level  trough  moved  onshore,  which  led  to  widespread 

precipitation  in  the  region. 

2 

February  09-10 

Quiescent  weather  was  in  place  with  a  500-hPa  ridge  centered 

over  central  California  at  12  UTC. 

3 

February  16-17 

An  upper-level  low  located  near  the  California/ Arizona 
border  with  Mexico  at  12  UTC  brought  precipitation  to  that 

portion  of  the  domain.  This  pattern  moved  south  and  east  over 
the  course  of  the  day. 

4 

March  01-02 

A  weak  shortwave  trough  resulted  in  precipitation  in  northern 
California  at  the  beginning  of  the  period  that  spread  to 

Nevada,  then  moved  southward  and  decreased  in  coverage. 

5 

March  05-06 

Widespread  high-level  cloudiness  due  to  weak  upper-level 
low  pressure  but  very  limited  precipitation. 

3.  Data  Preparation  using  MET  and  Visual  Basic 

We  used  MET  Point-Stat  to  compute  matched  pair  model  values  for  the  location 
of  each  surface  (2  and  10  meter)  observation  from  the  MADIS  dataset  in  each 
forecast  hour.  To  select  observations  closest  to  the  top  of  the  hour,  and  to 
eliminate  multiple  observations  for  each  hour,  we  set  the  duplicate  handling  of 
Point-Stat  to  SINGLE.  For  our  data,  we  extracted  the  following  fields:  site 
identification  code,  site  latitude,  site  longitude,  site  elevation,  the  measured  value 
of  the  meteorological  variable,  and  the  interpolated  forecast  value  from  the 
matched  pair  data  for  each  of  the  25  forecast  lead  times.  For  surface  observations, 
MET  derives  the  forecast  value  corresponding  to  a  given  observation  by 
horizontal  interpolation  of  the  surface  diagnostic  variable  (e.g.,  2-m  temperature, 
2-m  moisture,  10-m  winds)  from  the  enclosing  grid  points. 

MET  Point-Stat  produces  a  single  matched  pair  output  file  for  each  model  output 
hour  considered.  Although  this  file  is  useful  for  many  analyses,  it  is  far  easier  to 
combine  these  files  into  a  single  organized  structure  that  supports  GIS  analysis. 
We  organized  our  data  using  Microsoft  Access  to  compile  the  matched  pair  data 
for  each  forecast  hour  into  a  single  database,  and  we  used  custom  Visual  Basic  for 
Applications  scripts  to  organize  the  data  into  a  tabular  structure  that  can  be  easily 
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used  for  GIS  analysis.  Generically,  this  table  consists  of  4  columns  that  contain 
site  information,  a  column  containing  the  meteorological  variable  of  interest,  and 
then  a  sequence  of  columns  containing  the  forecast,  observed,  and  model  error 
(modeled-observed)  for  each  of  the  25  model  output  times. 

During  our  efforts  to  organize  our  data,  we  discovered  that  several  observation 
stations  from  different  locations  had  the  same  site  identification  code.  Because  the 
differences  were  as  much  as  a  kilometer  apart,  we  eliminated  these  station 
measurements  from  our  analysis.  We  further  restricted  our  analysis  to  the  stations 
that  reported  observations  in  each  of  the  25  h  of  a  given  case  day  (12  UTC  to  12 
UTC  the  next  day).  Our  intent  in  making  this  restriction  was  twofold:  1)  we  could 
make  hour-to-hour  comparisons  using  the  same  base  set  of  data,  and  2)  any 
remaining  stations  within  a  given  hour  would  serve  as  means  to  cross-validate  our 
analysis.  The  number  of  stations  considered  for  each  case  day  is  shown  in  Table 
2,  and  the  locations  of  these  stations  are  show  in  Figs.  2a-e,  where  a-e  correspond 
to  cases  1-5. 

Table  2  Count  of  stations  that  were  used  in  each  case  study  day,  for  measurements  made  2 
meters  and  10  meters  above  ground  level  (AGL) 


Altitude  (AGL) 


Case 

2  meters 

10  meters 

1 

135 

72 

2 

96 

48 

3 

122 

58 

4 

133 

80 

5 

69 

47 

The  MADIS  parameters  available  at  2  meters  AGL  include  temperature,  dew 
point,  and  relative  humidity.  The  available  10-meter  AGL  parameters  include  the 
“u”  and  “v”  components  of  the  wind,  as  well  as  the  magnitude  and  direction  of  the 
wind  speed.  For  this  study,  we  focused  our  attention  on  the  2-meter  observations 
of  relative  humidity  and  temperature.  We  will  consider  the  other  surface 
parameters — such  as  dew  point  and  wind — in  subsequent  studies. 
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Fig.  2a  The  location  of  the  weather  stations  that  contain  data  for  all  25  h  of  case  day  1 
overlaid  on  terrain  data  created  from  a  US  Geological  Survey  (USGS)  digital  elevation 
model  with  a  resolution  of  1/3  arc-second 
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Fig.  2b  The  location  of  the  weather  stations  that  contain  data  for  all  25  h  of  case  day  2 
overlaid  on  terrain  data  created  from  a  USGS  digital  elevation  model  with  a  resolution  of 
1/3  arc-second 
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Fig.  2c  The  location  of  the  weather  stations  that  contain  data  for  all  25  h  of  case  day  3 
overlaid  on  terrain  data  created  from  a  USGS  digital  elevation  model  with  a  resolution  of 
1/3  arc-second 
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Fig.  2d  The  location  of  the  weather  stations  that  contain  data  for  all  25  h  of  case  day  4 
overlaid  on  terrain  data  created  from  a  USGS  digital  elevation  model  with  a  resolution  of 
1/3  arc-second 
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Fig.  2e  The  location  of  the  weather  stations  that  contain  data  for  all  25  h  of  case  day  5 
overlaid  on  terrain  data  created  from  a  USGS  digital  elevation  model  with  a  resolution  of 
1/3  arc-second 

4.  Data  Analysis 


Our  first  step  was  to  examine  the  MET  model  performance  statistics  generated 
when  MET  created  the  matched-pair  output  data  for  the  inner  domain  (domain  3). 
In  Table  3,  we  summarize  these  statistics  for  2  h  in  each  of  the  5  cases;  4:00  PM 
Pacific  Standard  Time  (PST)  (0000  UTC)  to  capture  the  late  afternoon,  unstable 
Planetary  Boundary  Layer  (PBL)  conditions,  and  4:00  AM  PST  (1200  UTC)  on 
the  second  day  of  the  case  to  consider  the  early  morning,  stable  PBL  conditions.  It 
is  important  to  note  that  4:00  AM  is  also  the  last  forecast  produced  for  that  case 
day.  Appendix  A  expands  Table  3  to  include  other  hours  for  each  of  the  5  case 
days.  The  average  temperature  bias  error  in  Table  3  ranged  from  -3.55  to  2.42  K, 
and  the  average  relative  humidity  bias  error  ranged  from  -18.66  to  17.34  K, 
suggesting  that  the  model  perfonned  reasonably  well  over  domain  3.  However, 
because  MET  calculates  statistics  over  the  entire  domain,  large  positive  and 
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negative  bias  errors  could  be  cancelling  each  other  out  and  masking  poor  model 
perfonnance. 

Table  3  Mean  bias  (modeled-observed)  errors  for  each  of  the  5  case  days,  and  the  3 
meteorological  variables  modeled  at  the  Z2  surface  level  (2  meters  AGL):  4:00  PM  PST 
(0000  UTC)  and  4:00  AM  PST  (1200  UTC)  the  next  day.  These  statistics  were  calculated 
over  the  entire  domain. 


Case 

UTC 

Temperature 

(K) 

Dew  Point  (K) 

Relative 
Humidity  (%) 

1 

0000 

-0.17 

1.11 

5.83 

1200 

0.76 

1.40 

2.39 

2 

0000 

-2.89 

2.95 

12.27 

1200 

2.42 

-1.35 

-13.76 

3 

0000 

-0.65 

0.09 

1.94 

1200 

0.81 

-2.18 

-11.31 

4 

0000 

0.30 

0.23 

0.14 

1200 

1.58 

-2.18 

-18.66 

5 

0000 

-3.55 

3.48 

17.34 

1200 

-0.13 

-0.59 

-1.57 

Using  ArcGIS  version  10.2.2  (ESRI  2014),  we  created  quantile-quantile  (Q-Q) 
plots  of  the  model  errors  (modeled  value  minus  measured  value)  of  each  hour.  A 
Q-Q  plot  provides  a  visual  comparator  of  observed  or  experimental  data  to  a 
reference  probability  distribution  (Jolliffe  and  Stephenson  2011;  Wilks  2011).  We 
found  that  all  of  the  model  errors  were  approximately  normally  distributed; 
consequently,  we  used  conventional  statistical  techniques  that  are  predicated  on 
the  assumption  of  nonnality  in  our  analysis  of  the  model  errors. 

The  goal  of  this  study  was  to  develop  GIS  tools  for  detennining  whether  the 
terrain  influences  the  model  error.  Previous  studies  have  found  that  terrain 
variables  such  as  elevation  and  slope  have  a  strong  influence  on  climatological 
temperature  values  (Agnew  and  Palutikof  2000;  Brown  and  Comrie  2002; 
Chapman  and  Thornes  2003).  Other  candidate  explanatory  variables  include  land 
use  (Cheng  et  al.  2012)  and  soil  moisture  (Massey  et  al.  2014).  We  chose 
elevation  as  the  explanatory  terrain  variable  for  this  study  because  of  its  known 
affect  on  temperature  and  conducted  a  correlation  analysis  of  model  error  as  the 
dependent  variable  against  observation  station  elevation  as  the  independent 
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variable.  Future  studies  will  include  the  analysis  of  other  explanatory  variables 
such  as  slope,  aspect,  distance  from  the  coast,  and  land  use. 

4.1  Correlation  Statistics  for  the  Temperature  and  Relative 
Humidity  Bias  Errors  and  for  Terrain  Elevation 


We  performed  a  correlation  analysis  with  model  error  as  the  dependent  variable 
and  observation  elevation  as  the  independent  variable.  In  Figs.  3a-e,  where  a-e 
correspond  to  case  days  1-5,  we  plot  the  Pearson  correlation  coefficients  between 
the  relative  humidity  and  temperature  model  errors  as  dependent  variables  to 
elevation  as  the  independent  variable,  by  forecast  hour  for  each  case  study  day. 
We  tested  the  statistical  significance  of  each  correlation  coefficient  using  the 
Student  t-test  (Jolliffe  and  Stephenson  2011).  Correlation  coefficients  of  0.2  and 
greater  were  always  statistically  significant  with  the  number  of  measurements  that 
were  available.  On  case  study  days  with  a  larger  number  of  observations,  some 
correlation  coefficients  below  0.2  were  also  statistically  significant. 


2/7/12:  Correlation  of  Model  Bias  (Forecast-Observed)  with  Elevation 


Fig.  3a  The  Pearson  correlation  coefficient  computed  between  elevation  and  the  model 
errors  for  relative  humidity  (triangles)  and  surface  temperature  (circles)  for  each  hour  of 
case  day  1 
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2/9/12:  Correlation  of  Model  Bias  (Forecast-Observed)  with  Elevation 


Fig.  3b  The  Pearson  correlation  coefficient  computed  between  elevation  and  the  model 
errors  for  relative  humidity  (triangles)  and  surface  temperature  (circles)  for  each  hour  of 
case  day  2 
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2/16/12:  Correlation  of  Model  Bias  (Forecast-Observed)  with  Elevation 


Fig.  3c  The  Pearson  correlation  coefficient  computed  between  elevation  and  the  model 
errors  for  relative  humidity  (triangles)  and  surface  temperature  (circles)  for  each  hour  of 
case  day  3 
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3/1/12:  Correlation  of  Model  Bias  (Forecast-Observed)  with  Elevation 


Fig.  3d  The  Pearson  correlation  coefficient  computed  between  elevation  and  the  model 
errors  for  relative  humidity  (triangles)  and  surface  temperature  (circles)  for  each  hour  of 
case  day  4 
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3/5/12:  Correlation  of  Model  Bias  (Forecast-Observed)  with  Elevation 


Fig.  3e  The  Pearson  correlation  coefficient  computed  between  elevation  and  the  model 
errors  for  relative  humidity  (triangles)  and  surface  temperature  (circles)  for  each  hour  of 
case  day  5 

We  conducted  further  analysis  of  the  data  in  Fig.  3  to  discern  diurnal  trends  in  the 
correlation  coefficient  between  temperature  and  relative  humidity.  In  general, 
there  appear  to  be  no  strong  trends.  However,  we  noted  that  the  magnitude  of  the 
correlation  of  temperature  with  elevation  to  be  highest  in  the  afternoon  from 
12:00  PM  PST  to  4:00  PM  PST  except  on  February  7,  the  day  with  the  most 
widespread  precipitation.  For  relative  humidity  this  was  not  the  case.  The  hour  of 
highest  correlation  with  temperature  occurred  at  varying  times  depending  on  the 
case  study  day. 

Comparing  the  temporal  changes  in  the  correlation  coefficients  for  both  variables, 
we  noticed  that  there  appeared  to  be  an  inverse  relationship  between  those  of 
relative  humidity  and  temperature.  For  each  of  the  5  case  study  days,  we  found 
that  the  correlation  plots  of  temperature  and  relative  humidity  are  nearly  perfect 
mirror  images  of  each  other,  another  indication  of  the  strong  negative  correlation 
between  temperature  and  relative  humidity  forecast  errors.  Wherever  one  variable 
is  positively  correlated  with  elevation,  the  other  variable  is  negatively  correlated 
with  elevation. 


18 


This  obvious  pattern  raised  our  interest  in  determining  why  this  is  the  case.  We 
posed  the  hypothesis  that  the  relative  humidity  and  temperature  forecast  errors 
were  correlated  with  each  other  and  then  calculated  Pearson  correlation 
coefficients  between  temperature  error  and  relative  humidity  error  to  test  this 
hypothesis.  Figure  4  shows  the  plots  of  the  coefficients  by  forecast  hour  for  each 
case  study  day. 

With  the  exception  of  a  few  hours,  we  found  that  the  relative  humidity  and 
temperature  error  have  a  strong  negative  correlation  and  thus  are  not  independent. 
This  means  that  the  same  explanatory  variables  may  potentially  explain  the 
variance  in  the  forecast  error  of  both  parameters.  There  does  not  appear  to  be  a 
diurnal  pattern  to  the  variation  in  the  correlation  coefficient  between  temperature 
and  relative  humidity,  which  suggests  that  synoptic  conditions  could  partially 
account  for  some  of  the  variation  in  Fig.  4.  The  most  likely  explanation  for  the 
strong  inverse  correlation  between  temperature  and  relative  humidity  model  error 
lies  in  the  underlying  relationship  between  relative  humidity  and  temperature. 
Because  relative  humidity  is  a  function  of  the  temperature  and  water  vapor  mixing 
ratio  and  an  increasing  temperature  usually  means  decreasing  relative  humidity,  a 
negative  correlation  between  temperature  and  relative  humidity  error  is  to  be 
expected  even  if  the  terrain  is  flat. 


5AM  7  AM  9AM  11AM  1PM  3PM  5PM  7  PM  9PM  11PM  1AM  3AM 


Time  (Pacific  Standard  Time) 


Fig.  4  The  Pearson  correlation  coefficient  between  the  model  error  of  surface 
temperature  and  the  model  error  of  surface  relative  humidity 
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To  quantify  the  contribution  of  terrain  elevation  to  the  overall  variability  of  the 
model  error,  we  calculated  the  Pearson  correlation  coefficients  for  temperature 
and  relative  humidity  model  error  with  terrain  elevation.  The  percent  of  variance 
in  the  model  error  attributable  to  elevation  as  an  explanatory  variable  was  plotted 
by  output  hour  for  each  case  study  day.  Elevation  accounts  for  a  sometimes 
substantial  portion  of  the  variance  in  both  the  temperature  and  relative  humidity 
model  error  for  some  cases  and  times,  as  shown  in  Figs.  5a-e,  where  a-e  refer  to 
case  days  1-5.  On  average,  elevation  accounts  for  about  10%  to  20%  of  the  total 
variance  in  the  model  bias.  It  is  difficult  to  discern  a  diurnal  pattern,  but  elevation 
generally  accounts  for  the  greatest  percentage  of  the  model  bias  during  the 
afternoon  hours  from  12:00  PM  to  5:00  PM  PST.  For  example,  at  noon  on 
February  7,  2012  (case  day  1),  elevation  accounts  for  56%  of  the  temperature 
model  bias  and  72%  of  the  relative  humidity  model  bias  (elevation’s  strongest 
influence  among  the  5  case  days  investigated). 


2/7/2012:  Model  Bias  Variance  (Forecast-Observed)  Attributable  to  Elevation 


Hour  of  Day  (Pacific  Standard  Time) 


Fig.  5a  The  percentage  of  model  temperature  error  (circles)  and  in  the  relative  humidity 
error  (triangles)  variance  that  is  attributable  to  elevation  as  an  explanatory  variable  for  case 
day  1 
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2/9/2012:  Model  Bias  Variance  (Forecast-Observed)  Attributable  to  Elevation 
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Fig.  5b  The  percentage  of  model  temperature  error  (circles)  and  in  the  relative  humidity 
error  (triangles)  variance  that  is  attributable  to  elevation  as  an  explanatory  variable  for  case 
day  2 
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2/16/2012:  Model  Bias  Variance  (Forecast-Observed)  Attributable  to  Elevation 


Hour  of  Day  (Pacific  Standard  Time) 


Fig.  5c  The  percentage  of  model  temperature  error  (circles)  and  in  the  relative  humidity 
error  (triangles)  variance  that  is  attributable  to  elevation  as  an  explanatory  variable  for  case 
day  3 
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3/1/2012:  Model  Bias  Variance  (Forecast-Observed)  Attributable  to  Elevation 
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Fig.  5d  The  percentage  of  model  temperature  error  (circles)  and  in  the  relative  humidity 
error  (triangles)  variance  that  is  attributable  to  elevation  as  an  explanatory  variable  for  case 
day  4 
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Fig.  5e  The  percentage  of  model  temperature  error  (circles)  and  in  the  relative  humidity 
error  (triangles)  variance  that  is  attributable  to  elevation  as  an  explanatory  variable  for  case 
day  5 

For  each  case  day,  we  determined  the  number  of  model  output  hours  out  of  the  25 
available  (6  analysis  and  19  forecast  hours),  which  had  elevation  as  a  significant 
explanatory  variable.  Because  the  number  of  observation  stations  considered  for  a 
given  case  day  was  fixed  for  that  case  (see  Table  2),  we  could  determine  a 
threshold  (at  an  a  of  5%)  for  the  correlation  coefficient  above  which  the 
correlation  would  be  considered  significant.  Finally,  we  judged  that  if  a  large 
number  (more  than  half)  of  the  model  output  hours  had  significant  correlations, 
then  elevation  was  a  significant  explanatory  variable  for  that  variable  and  case. 
This  data  is  presented  in  Table  4  for  both  Temperature  and  Relative  Humidity. 
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Table  4  Statistics  related  to  the  percentage  of  variance  of  the  temperature  and  relative 
humidity  forecast  errors  explained  by  elevation 

Note:  For  both  threshold  and  maximum  R2,  the  values  are  the  Pearson  R2  coefficient 
multiplied  by  100  to  show  the  values  as  percentiles. 


Temperature 


Case 

Day 

Threshold  R2 
(%) 

No. 

Significant 

Output 

Hours 

Maximum  R2 
(%) 

Elevation  Significant? 

1 

4.19 

22 

55.71 

Yes 

2 

3.99 

14 

23.29 

Yes 

3 

2.79 

21 

40.18 

Yes 

4 

4.01 

18 

41.98 

Yes 

5 

5.43 

5 

20.18 

No 

Relative  Humidity 

Case 

Day 

Threshold  R2 
(%) 

No. 

Significant 

Output 

Hours 

Maximum  R2 
(%) 

Elevation  Significant? 

1 

2.95 

21 

72.30 

Yes 

2 

4.26 

9 

18.32 

No 

3 

2.84 

19 

38.68 

Yes 

4 

4.37 

19 

56.14 

Yes 

5 

5.71 

10 

23.29 

No 

4.2  Estimating  the  Forecast  Error  over  the  Entire  Domain  using 

GIS 

Using  a  GIS,  we  created  a  raster  mosaic  of  terrain  data  from  USGS  digital 
elevation  model  (DEM)  data  (USGS  2013).  The  USGS  DEM  has  a  resolution  of 
1/3  arc-second  and  a  datum  of  NAD  1983.  A  datum  is  a  mathematical  algorithm  to 
account  for  the  fact  that  the  earth  is  not  perfectly  spherical.  A  projection  corrects 
for  the  distortion  that  occurs  when  projecting  spherical  latitude  and  longitude 
coordinates  onto  the  flat  surface  of  paper  or  a  computer  screen.  WRF  uses  no 
datum  (assumes  that  the  earth  is  a  perfect  sphere),  and  we  used  the  Lambert  Conic 
Confonnal  projection  in  our  WRF  modeling. 

Matching  the  Point-Stat  matched  pair  data  with  the  DEM  proved  to  be  a  big 
challenge,  because  the  projection  and  datum  of  the  MADIS  data  were  not 
specified.  Upon  the  advice  of  Olga  Wilhehni  and  Jennifer  Boehnert  from  NCAR 
(2014),  we  used  the  NAD  1983  datum  without  a  projection  for  the  Point-Stat  data. 
With  these  GIS  settings,  the  Point-Stat  data  matched  the  USGS  DEM  very  well. 
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For  example,  weather  stations  located  on  the  tiny  Channel  Islands  off  the 
California  coast  were  correctly  geo-located  on  land  rather  than  the  Pacific  Ocean. 

To  estimate  the  model  errors  over  the  entire  domain,  we  treated  each  weather 
station  measurement  as  a  sample  point  from  a  random  field  and  employed  the 
Empirical  Bayesian  Kriging  (EBK)  tool  (Krivoruchko  2011;  2012)  of  ArcGIS  to 
estimate  a  continuous  error  surface.  These  surfaces  do  not  extend  to  the  western 
or  southern  extents  of  the  domain  because  there  is  insufficient  data  in  these  areas 
to  produce  a  statistically  significant  surface  estimate.  Contours  of  the  bias  error, 
corresponding  to  the  quintiles  of  the  error  distribution,  were  created  to  ensure  that 
the  underlying  terrain  was  visible. 

In  our  analysis  of  the  5  case  study  days,  we  studied  the  5  specific  hours  given  in 
Table  5.  We  chose  these  specific  times  to  examine  the  impact  of  the  data 
assimilation  and  the  stability  of  the  boundary  layer.  Under  clear  skies  and  weak 
synoptic  forcing,  at  5:00  AM  PST  (the  second  assimilation  hour)  we  would  expect 
the  boundary  layer  to  be  stable  because  the  sun  has  not  yet  warmed  the  surface  of 
the  earth.  Similarly  the  boundary  layer  is  often  unstable  at  about  0000  UTC  due  to 
a  continual  upward  heat  flux  following  the  earlier  maximum  solar  insolation 
period.  As  the  ground  cools  after  sunset,  the  boundary  layer  becomes  more  stable. 
Clear  sky  conditions  were  less  common  during  the  5  case  study  days  of  this  study. 
Precipitation  occurred  over  the  study  domain  on  case  days  1,  3,  and  4;  and  much 
of  the  domain  was  cloudy  on  case  day  5.  Under  cloudy  conditions  with 
precipitation,  we  expect  that  the  actual  tendency  would  have  been  for  more 
neutral/less  stable  conditions  during  nighttime  and  less  unstable  during  the  day. 

Table  5  Hours  chosen  for  further  GIS  study  and  the  rationale  for  these  choices 


Time  (h) 


Period 

UTC 

Forecast 

Lead 

Local 

(PST) 

Reason 

Analysis 

13 

18 

-5 

0 

5:00  AM 

10:00  AM 

First  hour  of  data  assimilation 
completed 

Beginning  of  forecast;  near  neutral 
boundary  layer 

0 

6 

4:00  PM 

Late  afternoon,  unstable  boundary  layer 

Forecast 

6 

12 

10:00  PM 

Early  evening,  slightly  stable  boundary 
layer 

12 

18 

4:00  AM 

Early  morning,  stable  boundary  layer 

The  effect  of  the  terrain  on  the  model  error  can  be  clearly  seen  in  the  March  1, 
2012  (case  day  4)  error  contour  maps  of  Fig.  6.  For  most  of  the  day,  the 
temperature  model  bias  is  positively  correlated  with  elevation,  while  relative 
humidity  model  bias  is  negatively  correlated  with  elevation.  However,  an 
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exception  is  5:00  AM  PST  (13  UTC;  Figure  6a,  top  row),  when  these  correlations 
are  reversed;  however,  the  correlation  coefficients  for  this  hour  are  not 
statistically  significant.  The  strong  inverse  relationship  between  temperature  and 
relative  humidity  is  clearly  evident  in  the  afternoon  hours  of  March  1  (Fig.  6b,  top 
row).  The  terrain  effects  can  also  be  seen  in  the  forecast  error  contour  maps  for 
February  9,  February  26,  and  March  5,  2012,  which  can  be  found  in  Appendix  B. 
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Fig.  6a  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for  March 
1,  2012  (case  day  4),  for  5:00  AM  PST  (13  UTC)  and  10:00  AM  PST  (18  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  5:00  AM  PST  (13  UTC,  1  h  into  preforecast)  and  the  bottom  row  is 
10:00  AM  PST  (18  UTC,  0-h  forecast). 
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Fig.  6b  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for  March 
1,  2012  (case  day  4),  for  4:00  PM  PST  (0  UTC)  and  10:00  PM  PST  (6  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  4:00  PM  PST  (0  UTC,  6-h  forecast)  and  the  bottom  row  is  10:00  PM 
PST  (6  UTC,  12-h  forecast). 
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Fig.  6c  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for  March 
1,  2012  (case  day  4),  for  4:00  AM  PST  (12  UTC)  on  March  2,  2012 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  This  row  is  for  4:00  AM  PST  the  next  day  (12  UTC,  18-h  forecast). 

We  also  created  error  contours  for  each  hour  of  February  7,  2012.  Selected  hours 
are  shown  in  Figs.  7a-g,  and  error  contours  for  every  hour  can  be  found  in 
Appendix  C.  From  10:00  AM  PST  to  1:00  PM  PST  (Figs.  7a  and  7b,  respectively) 
temperature  model  error  is  negatively  correlated  with  elevation,  while  the  relative 
humidity  model  error  is  positively  correlated.  These  correlations  are  evident  in 
both  the  error  contours  and  the  scatter  plots  of  forecast  error  versus  elevation. 
Between  5:00  PM  and  6:00  PM  PST  (Fig.  7c),  these  correlations  are  reversed  so 
that  temperature  forecast  error  is  positively  correlated  with  elevation,  while  the 
relative  humidity  forecast  error  is  negatively  correlated.  This  trend  continues  from 
7:00  PM  PST  to  10:00  PM  PST  (Figs.  7d  and  7e,  respectively).  Between  11:00 
PM  and  12:00  AM  PST  (Fig.  7f),  the  correlations  are  again  reversed. 
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Fig.  7a  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  10:00  AM  PST  (18  UTC)  and  11:00  AM  PST  (19  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  10:00  AM  PST  (18  UTC,  0-h  forecast)  and  the  bottom  row  is  11:00 
AM  PST  (19  UTC,  1-h  forecast). 
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Fig.  7b  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  12:00  AM  PST  (20  UTC)  and  1:00  AM  PST  (21  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  12:00  PM  PST  (20  UTC,  2-h  forecast)  and  the  bottom  row  is  1:00  PM 
PST  (21  UTC,  3-h  forecast). 
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Fig.  7c  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  5:00  PM  PST  (1  UTC)  and  6:00  PM  PST  (2  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  5:00  PM  PST  (1  UTC,  7-h  forecast)  and  the  bottom  row  is  6:00  PM 
PST  (2  UTC,  8-h  forecast). 
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Fig.  7d  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  7:00  PM  PST  (3  UTC)  and  8:00  PM  PST  (4  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  7:00  PM  PST  (3  UTC,  9-h  forecast)  and  the  bottom  row  is  8:00  PM 
PST  (4  UTC,  10-h  forecast). 
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Fig.  7e  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  9:00  PM  PST  (5  UTC)  and  10:00  PM  PST  (6  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  9:00  PM  PST  (5  UTC,  11-h  forecast)  and  the  bottom  row  is  10:00  PM 
PST  (6  UTC,  12-h  forecast). 
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Fig.  7f  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  11:00  PM  PST  (7  UTC)  and  12:00  AM  PST  (8  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  11:00  PM  PST  (7  UTC,  13-h  forecast)  and  the  bottom  row  is  12:00 
AM  PST  (8  UTC,  14-h  forecast). 
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Fig.  7g  EBK-created  isolines  of  constant  model  error  (forecast  minus  observed)  for 
February  7,  2012  (case  day  1)  for  1:00  AM  PST  (9  UTC)  and  2:00  AM  PST  (10  UTC) 

Note:  The  temperature  model  error  is  shown  in  the  left  column  and  the  relative  humidity 
model  error  in  the  right  column.  Each  map  shows  USGS  terrain  height  (shaded)  overlaid 
with  colored  isolines  showing  the  quintiles  of  the  error  distribution.  The  scatterplot  with 
each  map  shows  the  relationship  between  terrain  height  and  model  error  along  with  a  best 
fit  line.  The  top  row  is  1:00  AM  PST  (9  UTC,  15-h  forecast)  and  the  bottom  row  is  2:00  AM 
PST  (10  UTC,  16-h  forecast). 
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5.  Summary  and  Conclusions 


We  have  found  that  elevation  accounts  for  a  portion  of  the  variance  in  the 
temperature  and  relative  humidity  WRF  model  errors.  On  average,  elevation 
accounts  for  about  10%  to  20%  of  the  total  variance,  with  elevation  explaining 
more  than  50%  of  the  total  variance  during  the  afternoon  hours  of  February  7, 
2012.  By  contrast,  elevation  accounts  for  none  of  the  variance  during  the  early 
morning  hours  of  March  1  and  March  5.  Overall,  elevation  accounts  for  more 
variance  during  the  afternoon  hours  and  less  variance  during  the  early  morning 
hours;  but  we  cannot  distinguish  a  clear  diurnal  pattern  as  to  how  much  elevation 
affects  the  WRF  model  error. 

Future  studies  will  evaluate  the  effect  of  other  explanatory  terrain  variables  such 
as  slope,  aspect,  and  distance  from  the  coastline,  latitude,  longitude,  soil  moisture, 
and  land  use  on  bias  errors.  Further  research  is  needed  to  evaluate  why  terrain 
variables  such  as  elevation  affect  the  WRF  model  error.  It  would  also  be 
beneficial  to  evaluate  how  synoptic  weather  conditions  influence  the  relationship 
between  model  error  and  terrain  variables.  Analysis  of  weather  conditions  would 
likely  help  us  to  better  understand  the  day-to-day  (and  hour-to-hour)  variation  in 
the  model  bias  errors. 

Another  possible  source  of  model  error  may  be  geolocation  errors.  The  WRF 
model  assumes  that  the  earth  is  a  perfect  sphere  and  does  not  use  a  datum  to 
correct  for  the  fact  that  the  earth  is  not  perfectly  spherical.  Geolocation 
discrepancies  cause  errors  in  the  Coriolis  forcing  and  lateral  boundary  conditions 
as  well  as  in  the  locations  of  various  observations,  which  can  have  discrepancies 
of  as  much  as  20  km  depending  on  latitude;  a  condition  of  concern  for  local  area 
models  with  horizontal  resolutions  on  the  order  of  1  km  (Monaghan  et  al.  2013; 
see  especially  Fig.  1,  p.  2121). 

We  achieved  our  goal  for  this  study,  which  was  to  develop  analytical  methods  that 
go  beyond  the  types  of  conclusions  gained  from  domain-level,  aggregate 
statistics.  Using  the  tools  of  the  GIS,  the  model  errors  can  be  analyzed  to  show 
their  possible  dependence  on  terrain  characteristics  at  the  spatial  scales  of 
phenomena  of  interest  to  Army  Warfighters — such  as  mountain/valley  breezes 
and  land/sea  breezes.  This  study  demonstrates  that  GIS  tools  show  considerable 
promise  toward  our  goal  of  demonstrating  the  relationships  and  dependencies  of 
high-resolution  model  errors  to  terrain  characteristics  in  subdomains  where  the 
statistics  are  more  homogeneous. 
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Appendix  A.  Domain  Level  Errors  for  All  5  Case  Days 
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Table  A-l  Mean  bias  (modeled-observed)  errors  for  each  of  the  5  case  days,  and  the  3 
meteorological  variables  modeled  at  the  Z2  surface  level  (2  meters  AGL).  These  statistics 
were  calculated  over  the  entire  domain. 


March  5-6,2012 

Hour  (PST)  Hour(UTC)  Temperature  (K)  Dew  Point  Temp.  (K)  Relative  Humidity  (%) 


4:00  AM 

12 

0.81 

-0.33 

-2.51 

5:00  AM 

13 

-0.63 

-0.82 

-1.09 

10:00  AM 

18 

-0.03 

0.33 

1.31 

4:00  PM 

0 

-3.55 

3.48 

17.34 

10:00  PM 

6 

0.99 

-2.23 

-10.77 

4:00  AM 

12 

-0.13 

-0.59 

-1.57 
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Intentionally  Left  Blank. 
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Appendix  B.  Remaining  Case  Days  beyond  Fig.  6 


47 


Note:  Each  map  in  this  appendix  shows  US  Geological  Survey  (USGS)  terrain  height 
(shaded)  overlaid  with  colored  isolines  showing  the  quintiles  of  the  bias  distribution.  The  left 
panel  shows  a  scatterplot  of  the  relationship  between  terrain  height  and  model  bias  along 
with  a  best  fit  line. 


Fig.  B-la  Weather  Research  and  Forecasting  (WRF)  1-km  domain  model  surface  bias 
(forecast  minus  observed)  for  the  February  9,  2012,  case  for  temperature  (left  column)  and 
relative  humidity  (right  column).  The  top  row  is  5:00  AM  Pacific  Standard  Time  (PST)  (13 
coordinated  universal  time  [UTC],  1  h  into  preforecast)  and  the  bottom  row  is  10:00  AM 
PST  (18  UTC,  0-h  forecast). 
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Fig.  B-lb  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  February 
9,  2012,  case  for  temperature  (left  column)  and  relative  humidity  (right  column).  The  top 
row  is  4:00  PM  PST  (00  UTC,  6-h  forecast)  and  the  bottom  row  is  10:00  PM  PST  (06  UTC, 
12-h  forecast) 
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Fig.  B-lc  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  February 
9,  2012,  case  for  temperature  (left  figure)  and  relative  humidity  (right  figure),  4:00  AM  the 
next  day  (12  UTC,  18-h  forecast) 


50 


Porcenlile  Model 
Temperature  Bias  |K) 

i  10  AM  (IS  VTC) 

er>  ins 


Fig.  B-2a  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  February 
16,  2012,  case  for  temperature  (left  column)  and  relative  humidity  (right  column).  The  top 
row  is  5:00  AM  PST  (13  UTC,  1  h  into  preforecast)  and  the  bottom  row  is  10:00  AM  PST  (18 
UTC,  0-h  forecast). 
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Fig.  B-2b  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the 
February  16,  2012,  case  for  temperature  (left  column)  and  relative  humidity  (right  column). 
The  top  row  is  4:00  PM  PST  (00  UTC,  6-h  forecast)  and  the  bottom  row  is  10:00  PM  PST  (06 
UTC,  12-h  forecast). 
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Fig.  B-2c  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  February 
16,  2012,  case  for  temperature  (left  figure)  and  relative  humidity  (right  figure),  4:00  AM  the 
next  day  (12  UTC,  18-h  forecast) 
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Fig.  B-3a  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  March  5, 
2012,  case  for  temperature  (left  column)  and  relative  humidity  (right  column).  The  top  row 
is  5:00  AM  PST  (13  UTC,  1  h  into  preforecast)  and  the  bottom  row  is  10:00  AM  PST  (18 
UTC,  0-h  forecast). 
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Fig.  B-3b  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  March  5, 
2012,  case  for  temperature  (left  column)  and  relative  humidity  (right  column).  The  top  row 
is  4:00  PM  PST  (00  UTC,  6-h  forecast)  and  the  bottom  row  is  10:00  PM  PST  (06  UTC,  12-h 
forecast). 
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Fig.  B-3c  WRF  1-km  domain  model  surface  bias  (forecast  minus  observed)  for  the  March  5, 
2012,  case  for  temperature  (left  figure)  and  relative  humidity  (right  figure),  4:00  AM  the 
next  day  (12  UTC,  18-h  forecast) 
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Appendix  C.  Case  Day  1,  All  Hours 
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Fig.  C-l  February  7,  2012,  Weather  Research  and  Forecasting  (WRF)  error  (forecast  minus 
observed).  Elevation  versus  temperature  forecast  error  (left  column)  and  elevation  versus 
relative  humidity  forecast  error  (right  column).  The  top  row  is  4:00  AM  Pacific  Standard 
Time  (PST)  (12  UTC)  and  the  bottom  row  is  5:00  AM  PST  (13  UTC). 
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Fig.  C-2  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  6:00  AM  PST  (14  UTC),  and  the  bottom  row  is  7:00  AM 
PST  (15  UTC). 
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Fig.  C-3  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  8:00  AM  PST  (16  UTC),  and  the  bottom  row  is  9:00  AM 
PST  (17  UTC). 
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Fig.  C-4  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  10:00  AM  PST  (18  UTC),  and  the  bottom  row  is  11:00 
AM  PST  (19  UTC). 
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Fig.  C-5  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  12:00  PM  PST  (20  UTC),  and  the  bottom  row  is  1:00 
PM  PST  (21  UTC). 
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Fig.  C-6  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  2:00  PM  PST  (22  UTC),  and  the  bottom  row  is  3:00  PM 
PST  (23  UTC). 
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Fig.  C-7  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  4:00  PM  PST  (0  UTC),  and  the  bottom  row  is  5:00  PM 
PST  (1  UTC). 
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Fig.  C-8  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  6:00  PM  PST  (2  UTC),  and  the  bottom  row  is  7:00  PM 
PST  (3  UTC). 
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Fig.  C-9  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  8:00  PM  PST  (4  UTC),  and  the  bottom  row  is  9:00  PM 
PST  (5  UTC). 
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Fig.  C-10  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  10:00  PM  PST  (6  UTC),  and  the  bottom  row  is  11:00 
PM  PST  (7  UTC). 
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Fig.  C-ll  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  12:00  AM  PST  (8  UTC),  and  the  bottom  row  is  1:00  AM 
PST  (9  UTC). 


68 


Fig.  C-12  February  7,  2012,  WRF  error  (forecast  minus  observed).  Elevation  versus 
temperature  forecast  error  (left  column)  and  elevation  versus  relative  humidity  forecast 
error  (right  column).  The  top  row  is  2:00  AM  PST  (10  UTC),  and  the  bottom  row  is  3:00  AM 
PST  (11  UTC). 
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Fig.  C-13  February  7,  2012,  WRF  error  (forecast  minus  observed).  4:00  AM  (12  UTC)  the 
next  day,  elevation  versus  temperature  forecast  error  (left)  and  elevation  versus  relative 
humidity  forecast  error  (right). 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3-D 

ACARS 

AGL 

ARL 

DEM 

EBK 

FDDA 

GFS 

GIS 

h 

hPa 

K 

km 

LTC 

MAD  IS 

MET 

MODE 

MYJ 

MyWIDA 

NCAR 

NOAA 

NWP 

PBL 

PST 

Q-Q 


3 -Dimensional 

Aircraft  Communications  Addressing  and  Reporting  System 

above  ground  level 

US  Army  Research  Laboratory 

digital  elevation  model 

Empirical  Bayesian  Kriging 

Four  Dimensional  Data  Assimilation 

Global  Forecast  System 

Geographic  Information  System 

hour(s) 

hectopascal 

Degrees  Kelvin 

kilometer 

Lieutenant  Colonel 

Meteorological  Assimilation  Data  Ingest  System 

Model  Evaluation  Tools 

Method  for  Object-based  Diagnostic  Evaluation 

Mellor- Y  amada- J  anj  ic 

My  Weather  Impacts  Decision  Aid 

National  Center  for  Atmospheric  Research 

National  Oceanic  and  Atmospheric  Administration 

Numerical  Weather  Prediction 

Planetary  Boundary  Layer 

Pacific  Standard  Time 

quantile-quantile 
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RRTM 


Rapid  Radiative  Transfer  Model 


TAMDAR  Tropospheric  Airborne  Meteorological  Data  Reporting 

USGS  US  Geological  Survey 

UTC  coordinated  universal  time 

WRE-N  Weather  Running  Estimate-Nowcast 

WRF  Weather  Research  and  Forecasting 

WRF-ARW  Weather  Research  and  Forecasting-Advanced  Research  WRF 
WSMR  White  Sands  Missile  Range 
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